# Get Table 6
# =========================================================================================================
# remove everything
rm(list=ls())

# =================================================================================

# Load packages
library(foreign)

# Load the bootstrapped data
load("Boot_ss_Dt.RData")
load("Jack_ss_Dt.RData")

# ---------------------------------------------------------------------
# Get the bootstrapped CI.
B <- 400

# ================================================================================================
# ================================================================================================

# summary of the estimates
SPE <- summary(Bps1w[,1])
# summary of the Bootstrapped estimates
SPB <- apply(Bps1w, MARGIN = 2, summary)

# acceleration factor
SPE.jk <- apply(JKps1w, MARGIN = 2, summary)
acc.top <- apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# Q1: the 2nd position of SPE is the mean estimate
pos <- 2
dec <- 3
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q1 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
SPE[pos];ci.q1

# ---------------------------------------------------------------------
# Median: the 3rd position of SPE is the mean estimate
pos <- 3
dec <- 3
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
SPE[pos];ci.md

# ---------------------------------------------------------------------
# Q3: the 5th position of SPE is the mean estimate
pos <- 5
dec <- 3
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q3 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
SPE[pos];ci.q3

# ---------------------------------------------------------------------
# Percentage of positive SP obs.
dec <- 4
phio <- qnorm(apply((Bps1w[,1:400] < Bps1w[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKps1w, MARGIN = 1, mean) - JKps1w)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKps1w, MARGIN = 1, mean) - JKps1w)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(Bps1w)){
  temp <- as.numeric(quantile(Bps1w[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(Bps1w), dec), nsmall = dec)
per


# ================================================================================================
# ================================================================================================

# summary of the estimates
SPE <- summary(Bps2w[,1])
# summary of the Bootstrapped estimates
SPB <- apply(Bps2w, MARGIN = 2, summary)

# acceleration factor
SPE.jk <- apply(JKps2w, MARGIN = 2, summary)
acc.top <- apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# ---------------------------------------------------------------------
# Q1: the 2nd position of SPE is the mean estimate
pos <- 2
dec <- 3
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q1 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
SPE[pos];ci.q1

# ---------------------------------------------------------------------
# Median: the 3rd position of SPE is the mean estimate
pos <- 3
dec <- 3
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
SPE[pos];ci.md

# ---------------------------------------------------------------------
# Q3: the 5th position of SPE is the mean estimate
pos <- 5
dec <- 3
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q3 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
SPE[pos];ci.q3

# ---------------------------------------------------------------------
# Percentage of positive SP obs.
dec <- 4
phio <- qnorm(apply((Bps2w[,1:400] < Bps2w[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKps2w, MARGIN = 1, mean) - JKps2w)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKps2w, MARGIN = 1, mean) - JKps2w)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(Bps2w)){
  temp <- as.numeric(quantile(Bps2w[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(Bps2w), dec), nsmall = dec)
per

